Time: 0 seconds
# libraries
library(readr)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(gganimate)
library(gifski)
library(dplyr)
library(tidyr)
library(ggthemes)
library(ggpubr)
library(lme4)
library(ggeffects)
library(lmerTest)
library(car)
library(ggplot2)
library(survival)
library(coxme)
library(survminer)
library(ggfortify)
library(broom)
library(ehahelper)
library(dena)
library(tidyr)
Time: 0 seconds
sessionInfo()
R version 4.3.0 (2023-04-21)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.4
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: Europe/London
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] dena_0.1.0 ehahelper_0.3.9999 broom_1.0.4 ggfortify_0.4.16
[5] survminer_0.4.9 coxme_2.2-18.1 survival_3.5-5 car_3.1-2
[9] carData_3.0-5 ggeffects_1.2.2 ggpubr_0.6.0 ggthemes_4.2.4
[13] tidyr_1.3.0 dplyr_1.1.2 gifski_1.6.6-1 gganimate_1.0.8
[17] ggplot2_3.4.2 sjlabelled_1.2.0 sjmisc_2.8.9 sjPlot_2.8.14
[21] readr_2.1.4 bdsmatrix_1.3-6 lmerTest_3.1-3 lme4_1.1-33
[25] Matrix_1.5-4
loaded via a namespace (and not attached):
[1] gridExtra_2.3 rlang_1.1.1 magrittr_2.0.3 compiler_4.3.0
[5] mgcv_1.8-42 vctrs_0.6.2 stringr_1.5.0 pkgconfig_2.0.3
[9] crayon_1.5.2 fastmap_1.1.1 backports_1.4.1 labeling_0.4.2
[13] KMsurv_0.1-5 utf8_1.2.3 rmarkdown_2.21 tzdb_0.3.0
[17] haven_2.5.2 nloptr_2.0.3 purrr_1.0.1 bit_4.0.5
[21] xfun_0.39 cachem_1.0.8 jsonlite_1.8.4 progress_1.2.2
[25] tweenr_2.0.2 parallel_4.3.0 prettyunits_1.1.1 R6_2.5.1
[29] bslib_0.4.2 stringi_1.7.12 boot_1.3-28.1 jquerylib_0.1.4
[33] numDeriv_2016.8-1.1 estimability_1.4.1 Rcpp_1.0.10 knitr_1.42
[37] modelr_0.1.11 zoo_1.8-12 splines_4.3.0 tidyselect_1.2.0
[41] rstudioapi_0.14 abind_1.4-5 yaml_2.3.7 lattice_0.21-8
[45] tibble_3.2.1 withr_2.5.0 bayestestR_0.13.1 evaluate_0.20
[49] survMisc_0.5.6 pillar_1.9.0 insight_0.19.1 generics_0.1.3
[53] vroom_1.6.3 hms_1.1.3 munsell_0.5.0 scales_1.2.1
[57] minqa_1.2.5 xtable_1.8-4 glue_1.6.2 emmeans_1.8.5
[61] tools_4.3.0 data.table_1.14.8 ggsignif_0.6.4 forcats_1.0.0
[65] mvtnorm_1.1-3 cowplot_1.1.1 grid_4.3.0 colorspace_2.1-0
[69] nlme_3.1-162 eha_2.10.3 performance_0.10.3 cli_3.6.1
[73] km.ci_0.5-6 fansi_1.0.4 sjstats_0.18.2 gtable_0.3.3
[77] rstatix_0.7.2 sass_0.4.6 digest_0.6.31 farver_2.1.1
[81] htmltools_0.5.5 lifecycle_1.0.3 bit64_4.0.5 MASS_7.3-60
Time: 0 seconds
Data file ADD_I is available on the github page.
# read the data file "ADD_I"
ADD_I <- read_csv("Tables/ADD_I.csv", col_types = cols(Island = col_factor(levels = c("CN","AR", "CE", "DS", "FR"))))
# create the earliest born individuals on each island
Isfy <- ADD_I %>% select(founderyear,Island) %>% filter(!duplicated(Island))
# create a data frame with means and sd of each morphometric
ADD_mean <- ADD_I %>% group_by(Island,BirthYear) %>% summarise(RTsd = sd(RightTarsus,na.rm=TRUE),WLsd = sd(WingLength,na.rm=TRUE),BMsd = sd(BodyMass,na.rm=TRUE),HBsd=sd(HeadBill,na.rm=TRUE),RightTarsus=mean(RightTarsus,na.rm=TRUE),WingLength=mean(WingLength,na.rm=TRUE),BodyMass=mean(BodyMass,na.rm=TRUE),HeadBill=mean(HeadBill,na.rm=TRUE))
Time: 0 seconds
####RightTarsus ----
IRTlmer <- lmer(RightTarsus ~ Island*fYear + Sex + (1|Observer) + (1|Ageclass) + (1|BirdID), data = ADD_I)
summary(IRTlmer)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: RightTarsus ~ Island * fYear + Sex + (1 | Observer) + (1 | Ageclass) +
(1 | BirdID)
Data: ADD_I
REML criterion at convergence: 10511
Scaled residuals:
Min 1Q Median 3Q Max
-5.0560 -0.4445 0.0052 0.4324 5.5726
Random effects:
Groups Name Variance Std.Dev.
BirdID (Intercept) 0.3048654 0.55215
Observer (Intercept) 0.0562847 0.23724
Ageclass (Intercept) 0.0009875 0.03142
Residual 0.1407944 0.37523
Number of obs: 6154, groups: BirdID, 3258; Observer, 82; Ageclass, 5
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 2.417e+01 6.083e-02 2.713e+02 397.385 < 2e-16 ***
IslandAR -7.215e-02 1.686e-01 3.839e+03 -0.428 0.668783
IslandCE -2.510e-01 9.829e-02 3.925e+03 -2.553 0.010704 *
IslandDS 2.519e-01 8.080e-02 3.408e+03 3.117 0.001842 **
IslandFR 1.869e-01 1.236e-01 2.500e+03 1.512 0.130545
fYear 8.838e-03 2.530e-03 2.174e+03 3.493 0.000488 ***
SexMale 1.596e+00 2.244e-02 3.211e+03 71.130 < 2e-16 ***
IslandAR:fYear -6.760e-04 1.024e-02 3.539e+03 -0.066 0.947388
IslandCE:fYear 1.764e-03 6.196e-03 4.079e+03 0.285 0.775832
IslandDS:fYear -2.910e-03 5.840e-03 4.085e+03 -0.498 0.618312
IslandFR:fYear 1.554e-03 1.756e-02 1.981e+03 0.089 0.929486
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) IslnAR IslnCE IslnDS IslnFR fYear SexMal IsAR:Y IsCE:Y IsDS:Y
IslandAR -0.163
IslandCE -0.309 0.111
IslandDS -0.473 0.113 0.218
IslandFR -0.331 0.081 0.150 0.225
fYear -0.779 0.188 0.339 0.504 0.366
SexMale -0.177 -0.029 0.008 -0.033 -0.024 -0.016
IslndAR:fYr 0.149 -0.922 -0.096 -0.101 -0.076 -0.200 0.026
IslndCE:fYr 0.261 -0.094 -0.915 -0.173 -0.126 -0.323 -0.010 0.097
IslndDS:fYr 0.314 -0.065 -0.138 -0.849 -0.208 -0.379 0.008 0.069 0.123
IslndFR:fYr 0.102 -0.025 -0.050 -0.108 -0.850 -0.143 0.020 0.029 0.048 0.215
#predicting the model
IRTlmerdata <- ggpredict(IRTlmer,terms=c("fYear [all]","Island","Sex")) %>% rename(fYear=x,RightTarsus=predicted,Island=group,Sex=facet)
IRTlmerdata2 <- merge(IRTlmerdata,Isfy,by="Island", all.x=TRUE) %>% filter(!(Island== "AR" & fYear > 25)) %>% filter(!(Island== "CE" & fYear > 25)) %>% filter(!(Island== "DS" & fYear > 18)) %>% filter(!(Island== "FR" & fYear > 11))%>% mutate(BirthYear = fYear + founderyear) %>% group_by(Island,BirthYear) %>% summarise(RightTarsus = mean(RightTarsus, na.rm=TRUE), std.error = mean(std.error, na.rm=TRUE))
#plotting the model (shown as ggarrange below)
IRTmod <- ggplot(IRTlmerdata2, aes(x = BirthYear, y = RightTarsus, fill = Island)) +
geom_point(data = ADD_mean, mapping=aes(x=BirthYear,y=RightTarsus,colour=Island), size = 2) +
geom_errorbar(data=ADD_mean,aes(ymin=RightTarsus-RTsd,ymax=RightTarsus+RTsd, colour=Island),alpha=0.7) +
geom_smooth(aes(colour = Island),method = "loess", se = FALSE) +
geom_ribbon(aes(ymin=RightTarsus-std.error,ymax=RightTarsus+std.error),alpha=0.5) +
labs(x = "Birth Year", y= "Tarsus Length (mm)") +
xlim(1990,2022) +
theme_tufte(base_size = 15, base_family = "Arial") +
scale_color_colorblind() + scale_fill_colorblind()
Time: 0 seconds
#### Body Mass ----
IBMlmer <- lmer(BodyMass ~ Island*fYear + Sex + (1|Observer) + (1|Ageclass) + (1|BirdID), data = ADD_I)
summary(IBMlmer)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: BodyMass ~ Island * fYear + Sex + (1 | Observer) + (1 | Ageclass) +
(1 | BirdID)
Data: ADD_I
REML criterion at convergence: 16677.7
Scaled residuals:
Min 1Q Median 3Q Max
-3.6968 -0.5705 -0.0271 0.5440 4.2625
Random effects:
Groups Name Variance Std.Dev.
BirdID (Intercept) 0.24508 0.4951
Observer (Intercept) 0.05818 0.2412
Ageclass (Intercept) 0.14657 0.3828
Residual 0.53384 0.7306
Number of obs: 6554, groups: BirdID, 3398; Observer, 91; Ageclass, 5
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.438e+01 1.846e-01 5.339e+00 77.909 2.29e-09 ***
IslandAR 1.189e+00 1.570e-01 3.120e+03 7.577 4.63e-14 ***
IslandCE 2.792e-01 1.214e-01 3.666e+03 2.299 0.02158 *
IslandDS 7.385e-01 1.037e-01 1.096e+03 7.122 1.92e-12 ***
IslandFR 5.264e-01 1.668e-01 1.471e+03 3.156 0.00163 **
fYear -2.623e-03 3.107e-03 5.626e+02 -0.844 0.39893
SexMale 1.433e+00 2.635e-02 2.907e+03 54.405 < 2e-16 ***
IslandAR:fYear -7.870e-02 1.043e-02 2.346e+03 -7.546 6.38e-14 ***
IslandCE:fYear -3.708e-02 7.935e-03 4.479e+03 -4.673 3.05e-06 ***
IslandDS:fYear -5.069e-02 7.678e-03 2.097e+03 -6.601 5.14e-11 ***
IslandFR:fYear -5.400e-02 2.379e-02 1.097e+03 -2.270 0.02343 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) IslnAR IslnCE IslnDS IslnFR fYear SexMal IsAR:Y IsCE:Y IsDS:Y
IslandAR -0.071
IslandCE -0.107 0.098
IslandDS -0.176 0.107 0.184
IslandFR -0.113 0.071 0.120 0.187
fYear -0.307 0.188 0.294 0.462 0.320
SexMale -0.068 -0.016 0.007 -0.035 -0.027 -0.017
IslndAR:fYr 0.065 -0.869 -0.082 -0.095 -0.067 -0.205 0.016
IslndCE:fYr 0.092 -0.086 -0.910 -0.144 -0.097 -0.283 -0.007 0.092
IslndDS:fYr 0.117 -0.064 -0.114 -0.847 -0.181 -0.345 0.015 0.069 0.098
IslndFR:fYr 0.035 -0.024 -0.039 -0.097 -0.868 -0.127 0.023 0.027 0.035 0.201
#predicting the model
IBMlmerdata <- ggpredict(IBMlmer,terms=c("fYear [all]","Island","Sex")) %>% rename(fYear=x,BodyMass=predicted,Island=group,Sex=facet)
IBMlmerdata2 <- merge(IBMlmerdata,Isfy,by="Island", all.x=TRUE) %>% filter(!(Island== "AR" & fYear > 25)) %>% filter(!(Island== "CE" & fYear > 25)) %>% filter(!(Island== "DS" & fYear > 18)) %>% filter(!(Island== "FR" & fYear > 11))%>% mutate(BirthYear = fYear + founderyear) %>% group_by(Island,BirthYear) %>% summarise(BodyMass = mean(BodyMass, na.rm=TRUE), std.error = mean(std.error, na.rm=TRUE))
#plotting the model
IBMmod <- ggplot(IBMlmerdata2, aes(x = BirthYear, y = BodyMass, fill = Island)) +
geom_point(data = ADD_mean, mapping=aes(x=BirthYear,y=BodyMass,colour=Island), size = 2) +
geom_errorbar(data=ADD_mean,aes(ymin=BodyMass-BMsd,ymax=BodyMass+BMsd, colour=Island),alpha=0.7) +
geom_smooth(aes(colour = Island),method = "lm", se = FALSE) +
geom_ribbon(aes(ymin=BodyMass-std.error,ymax=BodyMass+std.error),alpha=0.5) +
labs(x = "Birth Year", y= "Body Mass (g)") +
xlim(1990,2022) +
theme_tufte(base_size = 15, base_family = "Arial") +
scale_color_colorblind() + scale_fill_colorblind()
Time: 0 seconds
#### Wing Length ----
IWLlmer <- lmer(WingLength ~ Island*fYear + Sex + (1|Observer) + (1|Ageclass) + (1|BirdID), data = ADD_I)
summary(IWLlmer)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: WingLength ~ Island * fYear + Sex + (1 | Observer) + (1 | Ageclass) +
(1 | BirdID)
Data: ADD_I
REML criterion at convergence: 23050.6
Scaled residuals:
Min 1Q Median 3Q Max
-3.8213 -0.5350 0.0161 0.5335 4.0139
Random effects:
Groups Name Variance Std.Dev.
BirdID (Intercept) 0.95940 0.9795
Observer (Intercept) 0.40157 0.6337
Ageclass (Intercept) 0.08901 0.2983
Residual 1.26489 1.1247
Number of obs: 6510, groups: BirdID, 3408; Observer, 91; Ageclass, 5
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 6.592e+01 1.905e-01 1.616e+01 346.012 < 2e-16 ***
IslandAR 2.490e-01 2.795e-01 3.717e+03 0.891 0.37296
IslandCE -1.869e-01 2.098e-01 4.116e+03 -0.891 0.37310
IslandDS -5.162e-01 1.865e-01 2.916e+03 -2.767 0.00569 **
IslandFR 1.871e-01 2.876e-01 2.198e+03 0.650 0.51549
fYear -2.900e-02 5.796e-03 1.679e+03 -5.004 6.22e-07 ***
SexMale 2.589e+00 4.624e-02 3.176e+03 55.991 < 2e-16 ***
IslandAR:fYear -4.495e-02 1.871e-02 3.681e+03 -2.403 0.01632 *
IslandCE:fYear 2.268e-03 1.354e-02 4.598e+03 0.168 0.86696
IslandDS:fYear 4.375e-02 1.350e-02 4.156e+03 3.241 0.00120 **
IslandFR:fYear 2.893e-02 4.134e-02 2.172e+03 0.700 0.48407
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) IslnAR IslnCE IslnDS IslnFR fYear SexMal IsAR:Y IsCE:Y IsDS:Y
IslandAR -0.123
IslandCE -0.198 0.099
IslandDS -0.315 0.101 0.190
IslandFR -0.225 0.072 0.131 0.197
fYear -0.551 0.177 0.302 0.464 0.344
SexMale -0.118 -0.015 0.005 -0.033 -0.026 -0.015
IslndAR:fYr 0.115 -0.873 -0.083 -0.092 -0.070 -0.200 0.013
IslndCE:fYr 0.174 -0.088 -0.910 -0.152 -0.112 -0.299 -0.007 0.096
IslndDS:fYr 0.218 -0.061 -0.121 -0.851 -0.195 -0.359 0.012 0.068 0.107
IslndFR:fYr 0.074 -0.023 -0.043 -0.099 -0.857 -0.135 0.024 0.028 0.040 0.217
optimizer (nloptwrap) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.0760292 (tol = 0.002, component 1)
#predicting the model
IWLlmerdata <- ggpredict(IWLlmer,terms=c("fYear [all]","Island","Sex")) %>% rename(fYear=x,WingLength=predicted,Island=group,Sex=facet)
IWLlmerdata2 <- merge(IWLlmerdata,Isfy,by="Island", all.x=TRUE) %>% filter(!(Island== "AR" & fYear > 25)) %>% filter(!(Island== "CE" & fYear > 25)) %>% filter(!(Island== "DS" & fYear > 18)) %>% filter(!(Island== "FR" & fYear > 11)) %>% mutate(BirthYear = fYear + founderyear) %>% group_by(Island,BirthYear)%>% summarise(WingLength = mean(WingLength, na.rm=TRUE), std.error = mean(std.error, na.rm=TRUE))
#plotting the model
IWLmod <- ggplot(IWLlmerdata2, aes(x = BirthYear, y = WingLength, fill = Island)) +
geom_point(data = ADD_mean, mapping=aes(x=BirthYear,y=WingLength,colour=Island), size = 2) +
geom_errorbar(data=ADD_mean,aes(ymin=WingLength-WLsd,ymax=WingLength+WLsd, colour=Island),alpha=0.7) +
geom_smooth(aes(colour = Island),method = "lm", se = FALSE) +
geom_ribbon(aes(ymin=WingLength-std.error,ymax=WingLength+std.error),alpha=0.5) +
labs(x = "Birth Year", y= "Wing Length (mm)") +
xlim(1990,2022) +
theme_tufte(base_size = 15, base_family = "Arial") +
scale_color_colorblind() + scale_fill_colorblind()
Time: 0 seconds
#### Head + Bill ----
IHBlmer <- lmer(HeadBill ~ Island*fYear + Sex + (1|Observer) + (1|Ageclass) + (1|BirdID), data = ADD_I)
summary(IHBlmer)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: HeadBill ~ Island * fYear + Sex + (1 | Observer) + (1 | Ageclass) +
(1 | BirdID)
Data: ADD_I
REML criterion at convergence: 9601.4
Scaled residuals:
Min 1Q Median 3Q Max
-4.1841 -0.4234 0.0279 0.4642 5.4733
Random effects:
Groups Name Variance Std.Dev.
BirdID (Intercept) 0.22985 0.4794
Observer (Intercept) 0.04913 0.2217
Ageclass (Intercept) 0.04152 0.2038
Residual 0.14507 0.3809
Number of obs: 5937, groups: BirdID, 3120; Observer, 80; Ageclass, 5
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 3.588e+01 1.085e-01 7.815e+00 330.737 <2e-16 ***
IslandAR 4.153e-01 2.063e-01 1.204e+03 2.013 0.0443 *
IslandCE -1.041e-01 9.292e-02 3.741e+03 -1.121 0.2624
IslandDS 8.141e-02 7.652e-02 3.252e+03 1.064 0.2874
IslandFR -7.000e-02 1.160e-01 2.381e+03 -0.603 0.5464
fYear -1.237e-03 2.492e-03 2.573e+03 -0.496 0.6196
SexMale 8.420e-01 2.075e-02 3.016e+03 40.579 <2e-16 ***
IslandAR:fYear -2.977e-02 1.559e-02 5.975e+02 -1.910 0.0566 .
IslandCE:fYear -8.098e-03 5.841e-03 3.978e+03 -1.386 0.1657
IslandDS:fYear 1.404e-03 5.456e-03 3.918e+03 0.257 0.7969
IslandFR:fYear 1.196e-02 1.638e-02 1.925e+03 0.730 0.4655
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) IslnAR IslnCE IslnDS IslnFR fYear SexMal IsAR:Y IsCE:Y IsDS:Y
IslandAR -0.066
IslandCE -0.177 0.079
IslandDS -0.277 0.088 0.233
IslandFR -0.197 0.067 0.163 0.243
fYear -0.444 0.155 0.351 0.525 0.386
SexMale -0.089 -0.017 -0.001 -0.036 -0.029 -0.021
IslndAR:fYr 0.046 -0.944 -0.049 -0.062 -0.051 -0.133 0.009
IslndCE:fYr 0.152 -0.064 -0.918 -0.188 -0.138 -0.333 -0.002 0.047
IslndDS:fYr 0.187 -0.049 -0.147 -0.851 -0.220 -0.393 0.012 0.041 0.131
IslndFR:fYr 0.062 -0.021 -0.054 -0.113 -0.848 -0.148 0.023 0.019 0.050 0.218
#predicting the model
IHBlmerdata <- ggpredict(IHBlmer,terms=c("fYear [all]","Island","Sex")) %>% rename(fYear=x,HeadBill=predicted,Island=group,SexEstimate=facet)
IHBlmerdata2 <- merge(IHBlmerdata,Isfy,by="Island", all.x=TRUE) %>% filter(!(Island== "AR" & fYear > 22)) %>% filter(!(Island== "CE" & fYear > 25)) %>% filter(!(Island== "DS" & fYear > 18)) %>% filter(!(Island== "FR" & fYear > 11))%>% mutate(BirthYear = fYear + founderyear)%>% group_by(Island,BirthYear)%>% summarise(HeadBill = mean(HeadBill, na.rm=TRUE), std.error = mean(std.error, na.rm=TRUE))
#plotting the model
IHBmod <- ggplot(IHBlmerdata2, aes(x = BirthYear, y = HeadBill, fill = Island)) +
geom_point(data = ADD_mean, mapping=aes(x=BirthYear,y=HeadBill,colour=Island), size = 2) +
geom_errorbar(data=ADD_mean,aes(ymin=HeadBill-HBsd,ymax=HeadBill+HBsd, colour=Island),alpha=0.7) +
geom_smooth(aes(colour = Island),method = "lm", se = FALSE) +
geom_ribbon(aes(ymin=HeadBill-std.error,ymax=HeadBill+std.error),alpha=0.25) +
labs(x = "Birth Year", y= "Head + Bill Length (mm)") +
xlim(1990,2022) +
theme_tufte(base_size = 15, base_family = "Arial") +
scale_color_colorblind() + scale_fill_colorblind()
Time: 0 seconds
Imods <- ggarrange(IRTmod,IBMmod,IWLmod,IHBmod, common.legend=TRUE, labels = c("A","B","C","D"))
Imods
Time: 0 seconds
# make a non-redundant data set for lifespan and reproduction
ADD_Real <- ADD_I %>% filter(Island == "CN") %>% group_by(BirdID) %>% mutate(across(c(RightTarsus,BodyMass,WingLength,HeadBill),mean, na.rm=TRUE)) %>% filter(!duplicated(BirdID)) %>% select(RightTarsus,BodyMass,WingLength,HeadBill,BirdID,SexEstimate,BirthYear,Lifespan,SurvStatus,ReproductiveOutput) %>% mutate(SexEstimate = as.factor(SexEstimate))
# filter for missing values
ADD_tidy <- ADD_Real %>% filter(!is.na(RightTarsus) & !is.na(BodyMass) & !is.na(WingLength) & !is.na(HeadBill) & !is.na(SexEstimate) & !is.na(Lifespan) ) %>% mutate(BodyCondition = BodyMass/RightTarsus) %>% mutate(SexEstimate = as.factor(SexEstimate))
Time: 0 seconds
###Lifespan analysis with coxme----
lscox0 <- coxme(Surv(Lifespan,SurvStatus) ~ RightTarsus+ BodyMass + WingLength + HeadBill +SexEstimate + (1|BirthYear), data= ADD_tidy)
lscox1 <- coxme(Surv(Lifespan,SurvStatus) ~ RightTarsus+ poly(BodyMass,2) + WingLength + HeadBill +SexEstimate + (1|BirthYear), data= ADD_tidy)
lscox2 <- coxme(Surv(Lifespan,SurvStatus) ~ RightTarsus*poly(BodyMass,2) + WingLength + HeadBill +SexEstimate + (1|BirthYear), data= ADD_tidy)
lscox3 <- coxme(Surv(Lifespan,SurvStatus) ~ poly(BodyCondition,2) + WingLength + HeadBill +SexEstimate + (1|BirthYear), data= ADD_tidy)
lscox4 <- coxme(Surv(Lifespan,SurvStatus) ~ RightTarsus+ poly(BodyMass,2)*SexEstimate + WingLength + HeadBill + (1|BirthYear), data= ADD_tidy)
AIC(lscox0,lscox1,lscox2,lscox3,lscox4) # lscox4 best
BIC(lscox0,lscox1,lscox2,lscox3,lscox4) # lscox4 best
# lscox4 was selected based on the best AIC and BIC values
lscox4
Cox mixed-effects model fit by maximum likelihood
Data: ADD_tidy
events, n = 1366, 1366
Iterations= 10 75
NULL Integrated Fitted
Log-likelihood -8500.56 -8298.841 -8240.677
Chisq df p AIC BIC
Integrated loglik 403.44 9.00 0 385.44 338.46
Penalized loglik 519.77 35.45 0 448.86 263.81
Model: Surv(Lifespan, SurvStatus) ~ RightTarsus + poly(BodyMass, 2) * SexEstimate + WingLength + HeadBill + (1 | BirthYear)
Fixed coefficients
coef exp(coef) se(coef) z p
RightTarsus 0.01881809 1.018996e+00 0.04997409 0.38 7.1e-01
poly(BodyMass, 2)1 -10.05540461 4.295298e-05 2.88625418 -3.48 4.9e-04
poly(BodyMass, 2)2 5.35344959 2.113361e+02 1.84019707 2.91 3.6e-03
SexEstimate1 0.54043577 1.716755e+00 0.10887255 4.96 6.9e-07
WingLength 0.02294382 1.023209e+00 0.02277508 1.01 3.1e-01
HeadBill 0.00173566 1.001737e+00 0.05698445 0.03 9.8e-01
poly(BodyMass, 2)1:SexEstimate1 -13.89794735 9.208696e-07 3.49003613 -3.98 6.8e-05
poly(BodyMass, 2)2:SexEstimate1 2.33161533 1.029456e+01 2.56984891 0.91 3.6e-01
Random effects
Group Variable Std Dev Variance
BirthYear Intercept 0.7408611 0.5488751
Time: 0.01 seconds
####predict lifespan----
ADD_round <- ADD_Real %>%mutate(across(c(RightTarsus, BodyMass, HeadBill, WingLength), round, 0)) %>% distinct(RightTarsus, BodyMass, HeadBill, WingLength,SexEstimate) %>% filter(!is.na(RightTarsus) & !is.na(BodyMass) & !is.na(WingLength) & !is.na(HeadBill) & !is.na(SexEstimate))
ADD_end <- bind_rows(ADD_tidy,ADD_round)
rs <- predict_coxme(lscox4,newdata = ADD_end,type="risk", se.fit=TRUE)
rsend <- ADD_end %>% add_columns(rr = rs$fit[,1], se.fit = rs$se.fit[,1])
ADD_tail <- rsend %>% tail(nrow(ADD_round)) %>% mutate(SexEstimate = as.factor(SexEstimate)) %>% group_by(BodyMass,SexEstimate) %>% summarise(mean_rr=mean(rr),mean_se.fit=mean(se.fit))
# plotting coxme
coxmeplot <- ggplot(ADD_tail,aes(x=BodyMass,y=mean_rr)) +
geom_line(aes(colour=SexEstimate)) +
geom_ribbon(aes(fill=SexEstimate,ymin=mean_rr-mean_se.fit,ymax=mean_rr+mean_se.fit), alpha = 0.25) +
labs(x = "Body Mass (g)", y = "Mortality Risk Score") +
guides(fill=guide_legend(title="Sex"), colour=guide_legend(title="Sex")) +
theme_tufte(base_size = 15, base_family = "Arial") +
scale_color_colorblind(labels = c("Female", "Male")) +
scale_fill_colorblind(labels = c("Female", "Male"))
Time: 0 seconds
###Reproduction glm----
ROglmF <- glm(ReproductiveOutput ~ RightTarsus+ BodyMass + WingLength + HeadBill ,data = ADD_tidy, family = "poisson")
ROglmR <- glm(ReproductiveOutput ~ poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm0 <- glm(ReproductiveOutput ~ RightTarsus+ BodyMass + WingLength + HeadBill + Lifespan + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm1 <- glm(ReproductiveOutput ~ RightTarsus+ BodyMass + WingLength + HeadBill + poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm2 <- glm(ReproductiveOutput ~ poly(RightTarsus,2)+ BodyMass + WingLength + HeadBill + poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm3 <- glm(ReproductiveOutput ~ RightTarsus+ poly(BodyMass,2) + WingLength + HeadBill + poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm4 <- glm(ReproductiveOutput ~ RightTarsus+ BodyMass + poly(WingLength,2) + HeadBill + poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm5 <- glm(ReproductiveOutput ~ RightTarsus+ BodyMass + WingLength + poly(HeadBill,2) + poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
AIC(ROglmF,ROglmR,ROglm0,ROglm1,ROglm2,ROglm3,ROglm4,ROglm5) # ROglm4 slightly better than ROglm1 second
BIC(ROglmF,ROglmR,ROglm0,ROglm1,ROglm2,ROglm3,ROglm4,ROglm5) # ROglm1 best
#ROglm1 was selected based on AIC and BIC scores
vif(ROglm1)
GVIF Df GVIF^(1/(2*Df))
RightTarsus 3.252977 1 1.803601
BodyMass 2.553710 1 1.598033
WingLength 2.550319 1 1.596972
HeadBill 2.274822 1 1.508251
poly(Lifespan, 2) 1.212419 2 1.049333
SexEstimate 3.638709 1 1.907540
BirthYear 1.266770 1 1.125509
summary(ROglm1)
Call:
glm(formula = ReproductiveOutput ~ RightTarsus + BodyMass + WingLength +
HeadBill + poly(Lifespan, 2) + SexEstimate + BirthYear, family = "poisson",
data = ADD_tidy)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 55.943199 6.971746 8.024 1.02e-15 ***
RightTarsus -0.004860 0.032759 -0.148 0.88207
BodyMass 0.075016 0.028836 2.602 0.00928 **
WingLength 0.015426 0.015372 1.004 0.31561
HeadBill -0.002388 0.038924 -0.061 0.95108
poly(Lifespan, 2)1 23.108050 1.113788 20.747 < 2e-16 ***
poly(Lifespan, 2)2 -5.454150 0.702063 -7.769 7.93e-15 ***
SexEstimate1 -0.042798 0.066926 -0.639 0.52251
BirthYear -0.028461 0.003275 -8.692 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 2018.43 on 736 degrees of freedom
Residual deviance: 879.58 on 728 degrees of freedom
(629 observations deleted due to missingness)
AIC: 3146.9
Number of Fisher Scoring iterations: 5
Time: 0.01 seconds
####predict reproduction----
# grouping variables for glm plot
ROplotdata <- ADD_tidy %>% mutate(BodyMass = case_when(BodyMass > 16.5 ~ "16.5", BodyMass < 14.5 ~ "14.5",BodyMass < 16.5 & BodyMass > 14.5 ~ "15.5")) %>% filter(!is.na(BodyMass))
# predicting glm
ROglmdat <- ggpredict(ROglm1, terms = c("Lifespan","BodyMass")) %>% rename(Lifespan=x,ReproductiveOutput=predicted,BodyMass=group)
#plotting glm
ROglmplot <- ggplot(ROglmdat, aes(x=Lifespan,y=ReproductiveOutput)) +
geom_point(data = ROplotdata, mapping=aes(x=Lifespan,y=ReproductiveOutput,colour=BodyMass),position="jitter") +
geom_smooth(aes(colour=BodyMass), se = FALSE) +
geom_ribbon(aes(ymin=conf.low,ymax=conf.high,fill=BodyMass), alpha = 0.25) +
labs(x = "Lifespan (years)", y= "Total Number of Offspring") +
theme_tufte(base_size = 15, base_family = "Arial") +
scale_color_colorblind() + scale_fill_colorblind()
Time: 0 seconds
fitplot <- ggarrange(coxmeplot,ROglmplot, labels = c("A","B"))
fitplot
Time: 0 seconds
Time taken to run the whole script
Time difference of 37.255 secs
Time: 0 seconds
Thank you for going through the code and supporting open research. This paper would not have been possible without all authors contributing towards it. Massive thanks to everyone that has shared their excitement for this paper.